---
title: "\\Large Replication code"
author: ""
date: ""
abstract: ""
output: 
  pdf_document:
    fig_caption: true
    keep_tex: true
fontsize: 11pt
header-includes:
    - \usepackage{caption}
    - \usepackage{bbm}
---

```{r packages}
# packages----
library(bit64) # long integers
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
library(knitr) #knitting
library(pander) #pandering
library(rasterVis) # plotting rasters
library(GGally) # for correlation plots
library(geosphere) # spherical distance
library(Rcpp) # for C++
library(RcppArmadillo) # for C++ (and Conley errors)
library(matrixStats) # matrix calculations
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(viridis) # for colors
library(Cairo) # pdf device
library(gtable) # for gtable filter
library(devtools)
#devtools::install_github("wilkelab/cowplot")
library(cowplot)
library(captioner)
library(splines)
library(survey) # svycontrast
library(texreg)
library(RColorBrewer)

# set cores for parallel processing----
registerDoMC(10)
```

```{r data}
k <- readRDS('replication_data.rds')
```

```{r fig1, eval=FALSE}
# cumulative cases----
cases.plot <- ggplot(data=k[year(lockdown_date_UTC)!=2021 & cum_cases>=10 & date >= as.Date('2020/03/01')], aes(x=date, y=cum_cases, group=fips, color=as.numeric(first_case, origin=min(first_case)))) + 
  geom_line(size=0.25, alpha=1) +
  scale_y_log10(labels = function(x) format(x, scientific = FALSE)) +
  annotation_logticks(sides = "l") +
  scale_color_gradient2(low="#006FFB", mid="lightblue", high="#FB8C00", guide=FALSE, midpoint = as.numeric(quantile(k$first_case, na.rm=T, p=0.025))) +
  ylab("Cumulative U.S. COVID-19 Cases") +
  xlab("Date") +
  theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=20, angle=0, face="bold"),
          axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
          axis.text.y = element_text(size=13, angle=0, face="bold"),
          axis.text.x = element_text(size=13, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          axis.ticks.y = element_line(),
          legend.title=element_blank()
    )

# cumulative deaths----
deaths.plot <- ggplot(data=k[year(lockdown_date_UTC)!=2021 & cum_deaths>=5 & date >= as.Date('2020/03/01')], aes(x=date, y=cum_deaths, group=fips, color=as.numeric(first_death, origin=min(first_death)))) + 
  geom_line(size=0.25, alpha=1) +
  scale_y_log10(labels = function(x) format(x, scientific = FALSE)) +
  annotation_logticks(sides = "l") +
  scale_color_gradient2(low="#006FFB", mid="lightblue", high="#FB8C00", guide=FALSE, midpoint = as.numeric(quantile(k$first_death, na.rm=T, p=0.025))) +
  ylab("Cumulative U.S. COVID-19 Fatalities") +
  xlab("Date") +
  theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=20, angle=0, face="bold"),
          axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
          axis.text.y = element_text(size=13, angle=0, face="bold"),
          axis.text.x = element_text(size=13, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          axis.ticks.y = element_line(),
          legend.title=element_blank()
    )

# plot----
png(filename = 'figure1.png', width=15, height=8, units = "in", type="cairo", res=300)
grid.arrange(cases.plot, deaths.plot, ncol=2)
grid.text("a", x=unit(0.09, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.585, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
dev.off()
```

```{r fig2, eval=FALSE}
# stay-at-home histogram----
stay.home <- k[year(lockdown_date_UTC)!=2021, .(lockdown_date_UTC = unique(lockdown_date_UTC)), by=.(fips)]

stay.home.plot <- ggplot(data=stay.home, aes(x=lockdown_date_UTC)) + 
  geom_histogram(size=0.5, fill="#FB8C00", color="#006FFB", alpha=0.75, bins=22) +
  ylab("Number of Counties") +
  xlab("Date of Stay-at-Home Order") +
  theme(plot.title = element_text(face="bold", size=40),
          axis.title.x = element_text(size=20, angle=0, face="bold"),
          axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
          axis.text.y = element_text(size=13, angle=0, face="bold"),
          axis.text.x = element_text(size=13, face="bold"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          axis.line.x = element_line(),
          axis.line.y = element_line(),
          axis.ticks.y = element_line(),
          legend.title=element_blank()
    )
  
# lockdown means plot----
lockdown.means <- k[abs(days_since_lockdown) <= (26), .(cases_growth = frollmean(cases_growth, 7, na.rm = T)), by = days_since_lockdown]
lockdown.means <- lockdown.means[, .(cases_growth = mean(cases_growth, na.rm=T)), by=.(days_since_lockdown)]

lockdown.plot <- ggplot(data=lockdown.means, aes(x=days_since_lockdown, y=cases_growth)) + 
  geom_point(size=4, fill="#FB8C00", color="white", pch=21, alpha=1) +
  geom_vline(xintercept=0, size=1, linetype="dotted", color="gray80") +
  ylab("Weekly Growth Rate of Cases (Percentage Points)") +
  xlab("Days Since Stay-at-Home") +
  coord_cartesian(ylim=c(-1,18), expand=TRUE) +
  theme(plot.title = element_text(face="bold", size=40),
        axis.title.x = element_text(size=20, angle=0, face="bold"),
        axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
        axis.text.y = element_text(size=13, angle=0, face="bold"),
        axis.text.x = element_text(size=13, face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line(),
        legend.title=element_blank()
  )

# plot----
png(filename = 'figure2.png', width=15, height=8, units = "in", type="cairo", res=300)
grid.arrange(stay.home.plot, lockdown.plot, ncol=2)
grid.text("a", x=unit(0.07, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.565, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
dev.off()
```

```{r fig3, eval=FALSE}
# create data for plotting----

# graph of cumulative cases by day by date of lockdown.
k[,the.date := as.Date(date)]
k[,the.lock.date := as.Date(lockdown_date_UTC)]
k[,the.reopen := as.Date(reopen_date_UTC)]
if (T) {
  # Collapse lockdown dates into pairs to smooth plotting.
  locks = k[,data.table(the.lock.date=sort(unique(the.lock.date)))]
  locks[,the.lock.2 := data.table::shift(the.lock.date,-1)]
  locks[,the.lock.3 := ifelse(seq_len(.N)%%2==0,paste(the.lock.date),paste(the.lock.2))]
  locks[,the.lock.2 := NULL]
  k = merge(k,locks,by="the.lock.date",all.x=T)
  k[,the.lock.date := as.Date(the.lock.3)]
}

# aggregate to the date-lockdown.type, include post-reopen.
AGG = k[,lapply(.SD,sum),keyby=c("the.lock.date","the.date"),.SDcols=c("cases","cum_cases","deaths","cum_deaths")] # tests only at state level ,"tests","cum_tests")]

# calculate lags of log cumulative cases.
AGG[cum_cases > 0,l.c.cases := log(cum_cases)]
setkey(AGG,the.lock.date,the.date)
for (l in c(1,2,3,4)) {
  AGG[,sprintf("L%s.l.c.cases",l) := data.table::shift(l.c.cases,l),by="the.lock.date"]
}

# rolling sum of new cases
AGG[,weekly.cases := 7*frollmean(cases,7),by="the.lock.date"]

# rolling sum of new deaths
AGG[,weekly.deaths := 7*frollmean(deaths,7),by="the.lock.date"]

# segments for stay-at-home dates
seg.cases <- AGG[the.lock.date == the.date, .(x = the.lock.date, 
                                        xend = the.lock.date, 
                                        y = 0, 
                                        yend = weekly.cases), by=the.lock.date]

seg.deaths <- AGG[the.lock.date == the.date, .(x = the.lock.date, 
                                        xend = the.lock.date, 
                                        y = 0, 
                                        yend = weekly.deaths), by=the.lock.date]

# plots----
weeklycases <- ggplot() + 
            geom_segment(data=seg.cases, aes(x=x, 
                                       xend=xend, 
                                       y=y,
                                       yend=yend,
                                       group=the.lock.date),
                         linetype="dotted", color="gray50") +
            geom_line(data=AGG[cum_cases >= 50 & !is.na(the.lock.date),], 
                      aes(x=the.date, y=weekly.cases, group=the.lock.date, color="lightblue"),
                      size=1) +
            geom_line(data=AGG[cum_cases >= 50 & is.na(the.lock.date),], 
                      aes(x=the.date, y=weekly.cases, group=the.lock.date, color="#FB8C00"),
                      size=2) +
            geom_point(data=seg.cases, aes(x=x, y=yend), color="gray50", size=1) +
            scale_y_log10(labels = function(x) format(x, scientific = FALSE)) +
            annotation_logticks(sides = "l") +
            scale_color_manual(values = c("lightblue","#FB8C00"),
                               breaks = c("lightblue","#FB8C00"),
                               labels = c("Stay-at-Home","No Stay-at-Home")) +
            annotate("text", x = as.Date('2020/04/18'), y = 200, label = "\u2190 Date Stay-at-Home\nEnacted", size=5) +
            ylab("Weekly COVID-19 Cases") +
            xlab("Date") +
            theme(plot.title = element_text(face="bold", size=40),
                    axis.title.x = element_text(size=20, angle=0, face="bold"),
                    axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
                    axis.text.y = element_text(size=13, angle=0, face="bold"),
                    axis.text.x = element_text(size=13, face="bold"),
                    panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    axis.line.x = element_line(),
                    axis.line.y = element_line(),
                    axis.ticks.y = element_line(),
                    legend.title=element_blank(),
                    legend.position = c(.2, 0.9),
                    legend.text=element_text(size=14, face="bold"),
                    legend.key = element_rect(colour = "black", fill="transparent")
              )

weeklydeaths <- ggplot() + 
            geom_segment(data=seg.deaths, aes(x=x, 
                                       xend=xend, 
                                       y=y,
                                       yend=yend,
                                       group=the.lock.date),
                         linetype="dotted", color="gray50") +
            geom_line(data=AGG[cum_deaths >= 5 & !is.na(the.lock.date),], 
                      aes(x=the.date, y=weekly.deaths, group=the.lock.date, color="lightblue"),
                      size=1) +
            geom_line(data=AGG[cum_deaths >= 5 & is.na(the.lock.date),], 
                      aes(x=the.date, y=weekly.deaths, group=the.lock.date, color="#FB8C00"),
                      size=2) +
            geom_point(data=seg.deaths, aes(x=x, y=yend), color="gray50", size=1) +
            scale_y_log10(labels = function(x) format(x, scientific = FALSE)) +
            annotation_logticks(sides = "l") +
            scale_color_manual(values = c("lightblue","#FB8C00"),
                               breaks = c("lightblue","#FB8C00"),
                               labels = c("Stay-at-Home","No Stay-at-Home")) +
            annotate("text", x = as.Date('2020/04/18'), y = 16, label = "\u2190 Date Stay-at-Home\nEnacted", size=5) +
            ylab("Weekly COVID-19 Fatalities") +
            xlab("Date") +
            theme(plot.title = element_text(face="bold", size=40),
                    axis.title.x = element_text(size=20, angle=0, face="bold"),
                    axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
                    axis.text.y = element_text(size=13, angle=0, face="bold"),
                    axis.text.x = element_text(size=13, face="bold"),
                    panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    axis.line.x = element_line(),
                    axis.line.y = element_line(),
                    axis.ticks.y = element_line(),
                    legend.title=element_blank(),
                    legend.position = c(.2, 0.9),
                    legend.text=element_text(size=14, face="bold"),
                    legend.key = element_rect(colour = "black", fill="transparent")
              )

# plot----
png(filename = 'figure3.png', width=15, height=8, units = "in", type="cairo", res=300)
grid.arrange(weeklycases, weeklydeaths, ncol=2)
grid.text("a", x=unit(0.095, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.59, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
dev.off()
```

```{r fig4, eval=FALSE}
# get data and results for plotting----
k <- readRDS("replication_data.rds")
k[is.na(cum_tests), cum_tests := 0]
p <- data.table(read.csv("co-est2019-alldata.csv"))
p[, fips := as.character(STATE * 1000 + COUNTY)]
k <- merge(k, p[, .(fips, POPESTIMATE2019)], all.x = T, by = "fips")

results <- results2 <- NULL
for(i in (-14) : 26) {
    print(i)
    d <- k[, .(cum_cases = sum(cum_cases), 
               cum_tests = mean(cum_tests), 
               cum_deaths = sum(cum_deaths), 
               weight = .N), 
           by = .(date = as.Date(date), 
                  lockdown_date = as.Date(lockdown_date_UTC))]
    n <- k[is.na(lockdown_date_UTC), 
           .(cum_cases = sum(cum_cases), 
             cum_tests = mean(cum_tests), 
             cum_deaths = sum(cum_deaths), 
             weight = .N), 
           by = .(date = as.Date(date))]
    d0 <- d[lockdown_date == date][order(lockdown_date)]
    d07 <- d[lockdown_date - 7 == date][order(lockdown_date)]
    d014 <- d[lockdown_date - 14 == date][order(lockdown_date)]
    d00 <- d[lockdown_date + i - 14 == date][order(lockdown_date)]
    d7 <- d[lockdown_date + i - 7 == date][order(lockdown_date)]
    d14 <- d[lockdown_date + i == date][order(lockdown_date)]

    d0[, cum_cases_07 := d07$cum_cases]
    d0[, cum_cases_014 := d014$cum_cases]
    d0[, cum_cases_0 := d00$cum_cases]
    d0[, cum_cases_7 := d7$cum_cases]
    d0[, cum_cases_14 := d14$cum_cases]

    d0[, cum_tests_07 := d07$cum_tests]
    d0[, cum_tests_014 := d014$cum_tests]
    d0[, cum_tests_0 := d00$cum_tests]
    d0[, cum_tests_7 := d7$cum_tests]
    d0[, cum_tests_14 := d14$cum_tests]

    d0[, cum_deaths_07 := d07$cum_deaths]
    d0[, cum_deaths_014 := d014$cum_deaths]
    d0[, cum_deaths_0 := d00$cum_deaths]
    d0[, cum_deaths_7 := d7$cum_deaths]
    d0[, cum_deaths_14 := d14$cum_deaths]

    d0[, n_weight := n[match(d0$date, n$date), weight]]

    d0[, n_cum_cases := n[match(d0$date,n$date), cum_cases]]
    d0[, n_cum_cases_07 := n[match(d0$date - 7, n$date), cum_cases]]
    d0[, n_cum_cases_014 := n[match(d0$date - 14, n$date), cum_cases]]
    d0[, n_cum_cases_0 := n[match(d0$date + i - 14, n$date), cum_cases]]
    d0[, n_cum_cases_7 := n[match(d0$date + i - 7, n$date), cum_cases]]
    d0[, n_cum_cases_14 := n[match(d0$date + i, n$date), cum_cases]]

    d0[, n_cum_tests := n[match(d0$date,n$date), cum_tests]]
    d0[, n_cum_tests_07 := n[match(d0$date - 7, n$date), cum_tests]]
    d0[, n_cum_tests_014 := n[match(d0$date - 14, n$date), cum_tests]]
    d0[, n_cum_tests_0 := n[match(d0$date + i - 14, n$date), cum_tests]]
    d0[, n_cum_tests_7 := n[match(d0$date + i - 7, n$date), cum_tests]]
    d0[, n_cum_tests_14 := n[match(d0$date + i, n$date), cum_tests]]

    d0[, n_cum_deaths := n[match(d0$date,n$date), cum_deaths]]
    d0[, n_cum_deaths_07 := n[match(d0$date - 7, n$date), cum_deaths]]
    d0[, n_cum_deaths_014 := n[match(d0$date - 14, n$date), cum_deaths]]
    d0[, n_cum_deaths_0 := n[match(d0$date + i - 14, n$date), cum_deaths]]
    d0[, n_cum_deaths_7 := n[match(d0$date + i - 7, n$date), cum_deaths]]
    d0[, n_cum_deaths_14 := n[match(d0$date + i, n$date), cum_deaths]]

    d0[, prior_change := log(cum_cases - cum_cases_07 + 1) - log(cum_cases_07 - cum_cases_014 + 1)]
    d0[, posterior_change := log(cum_cases_14 - cum_cases_7 + 1) - log(cum_cases_7 - cum_cases_0 + 1)]
    d0[, n_prior_change := log(n_cum_cases - n_cum_cases_07 + 1) - log(n_cum_cases_07 - n_cum_cases_014 + 1)]
    d0[, n_posterior_change := log(n_cum_cases_14 - n_cum_cases_7 + 1) - log(n_cum_cases_7 - n_cum_cases_0 + 1)]
    d0[, change_in_change := posterior_change - prior_change]
    d0[, n_change_in_change := n_posterior_change - n_prior_change]

    d0[, prior_change_tests := log(cum_tests - cum_tests_07 + 1) - log(cum_tests_07 - cum_tests_014 + 1)]
    d0[, posterior_change_tests := log(cum_tests_14 - cum_tests_7 + 1) - log(cum_tests_7 - cum_tests_0 + 1)]
    d0[, n_prior_change_tests := log(n_cum_tests - n_cum_tests_07 + 1) - log(n_cum_tests_07 - n_cum_tests_014 + 1)]
    d0[, n_posterior_change_tests := log(n_cum_tests_14 - n_cum_tests_7 + 1) - log(n_cum_tests_7 - n_cum_tests_0 + 1)]
    d0[, change_in_change_tests := posterior_change_tests - prior_change_tests]
    d0[, n_change_in_change_tests := n_posterior_change_tests - n_prior_change_tests]

    d0[, prior_change_deaths := log(cum_deaths - cum_deaths_07 + 1) - log(cum_deaths_07 - cum_deaths_014 + 1)]
    d0[, posterior_change_deaths := log(cum_deaths_14 - cum_deaths_7 + 1) - log(cum_deaths_7 - cum_deaths_0 + 1)]
    d0[, n_prior_change_deaths := log(n_cum_deaths - n_cum_deaths_07 + 1) - log(n_cum_deaths_07 - n_cum_deaths_014 + 1)]
    d0[, n_posterior_change_deaths := log(n_cum_deaths_14 - n_cum_deaths_7 + 1) - log(n_cum_deaths_7 - n_cum_deaths_0 + 1)]
    d0[, change_in_change_deaths := posterior_change_deaths - prior_change_deaths]
    d0[, n_change_in_change_deaths := n_posterior_change_deaths - n_prior_change_deaths]

    r <- data.table(
    log_weekly_change = c(d0[, posterior_change], d0[, prior_change], d0[, n_posterior_change], d0[, n_prior_change]),
    log_weekly_change_tests = c(d0[, posterior_change_tests], d0[, prior_change_tests], d0[, n_posterior_change_tests], d0[, n_prior_change_tests]),
    log_weekly_change_deaths = c(d0[, posterior_change_deaths], d0[, prior_change_deaths], d0[, n_posterior_change_deaths], d0[, n_prior_change_deaths]),
    county = c(d0[, lockdown_date], d0[, lockdown_date], rep(as.Date("2021-01-01"), dim(d0)[1]), rep(as.Date("2021-01-01"), dim(d0)[1])),
    date = rep(d0[, lockdown_date], 4),
    after_lockdown = c( rep(1, dim(d0)[1]), rep(0, dim(d0)[1]), rep(1, dim(d0)[1]), rep(0, dim(d0)[1]) ),
    is_lockdown_county = c( rep(1, dim(d0)[1]), rep(1, dim(d0)[1]), rep(0, dim(d0)[1]), rep(0, dim(d0)[1]) ),
    weight = c(d0[, weight], d0[, weight], d0[, n_weight], d0[, n_weight])
    )
    
    # drop county where weekly results in one county briefly go backwards due to revision
    r <- r[!is.infinite(log_weekly_change)]
    
    out <- summary(felm(log_weekly_change ~ log_weekly_change_tests + after_lockdown * is_lockdown_county | date | 0 | county, weights = r$weight, data = r, psdef=F))
    results <- rbind(results, c(i,out$coefficients[dim(out$coefficients)[1],]))

    out2 <- summary(felm(log_weekly_change_deaths ~ log_weekly_change_tests + after_lockdown * is_lockdown_county | date | 0 | county, weights = r$weight, data = r, psdef=F))
    results2 <- rbind(results2, c(i,out2$coefficients[dim(out2$coefficients)[1],]))
}

# set up data for plotting----
cases <- as.data.table(results)
setnames(cases, c('stay_day','coef','se','t','p'))
cases <- cases[, .(stay_day, coef, se)]
cases[stay_day==0, ":=" (coef=0, se=0)]
cases[, hi := coef + 1.96*se]
cases[, lo := coef - 1.96*se]

deaths <- as.data.table(results2)
setnames(deaths, c('stay_day','coef','se','t','p'))
deaths <- deaths[, .(stay_day, coef, se)]
deaths[stay_day==0, ":=" (coef=0, se=0)]
deaths[, hi := coef + 1.96*se]
deaths[, lo := coef - 1.96*se]

# felm-based plot----
casesplot <- ggplot() +
  geom_hline(aes(yintercept=0), size=1, linetype="dotted", color='gray80') +
  geom_vline(aes(xintercept=0), size=1, linetype="dotted", color='gray80') +
  geom_ribbon(data=cases, aes(x=stay_day, ymin=lo, ymax=hi), fill="transparent", alpha=0.1) +
  geom_line(data=cases, aes(x=stay_day, y=lo), linetype="dotted", color="#006FFB", size=1, alpha=1) +
  geom_line(data=cases, aes(x=stay_day, y=hi), linetype="dotted", color="#006FFB", size=1, alpha=1) +
  geom_line(data=cases, aes(x=stay_day, y=coef), color="#FB8C00", size=1, alpha=1) +
  geom_point(data=cases, aes(x=stay_day, y=coef), size=4, fill="#FB8C00", color="white", pch=21, alpha=1) +
  ylab("\u0394 Log Weekly Cases\nin Counties With Stay-at-Home vs. Counties Without") +
  xlab("Days Since Stay-at-Home") +
  coord_cartesian(xlim=c(min(cases$stay_day), max(cases$stay_day)), ylim=c(min(deaths$lo), max(deaths$hi)), expand=TRUE) +
  theme(
        axis.title.x = element_text(size=20, angle=0, face="bold"),
        axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
        axis.text.y = element_text(size=13, angle=0, face="bold"),
        axis.text.x = element_text(size=13, face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line(),
        legend.title=element_blank(),
        legend.position = c(.2, 0.9),
        legend.text=element_text(size=14, face="bold"),
        legend.key = element_rect(colour = "black", fill="transparent")
)

deathsplot <- ggplot() +
  geom_hline(aes(yintercept=0), size=1, linetype="dotted", color='gray80') +
  geom_vline(aes(xintercept=0), size=1, linetype="dotted", color='gray80') +
  geom_ribbon(data=deaths, aes(x=stay_day, ymin=lo, ymax=hi), fill="transparent", alpha=0.1) +
  geom_line(data=deaths, aes(x=stay_day, y=lo), linetype="dotted", color="#006FFB", size=1, alpha=1) +
  geom_line(data=deaths, aes(x=stay_day, y=hi), linetype="dotted", color="#006FFB", size=1, alpha=1) +
  geom_line(data=deaths, aes(x=stay_day, y=coef), color="#FB8C00", size=1, alpha=1) +

  geom_point(data=deaths, aes(x=stay_day, y=coef), size=4, fill="#FB8C00", color="white", pch=21, alpha=1) +
  ylab("\u0394 Log Weekly Fatalities\nin Counties With Stay-at-Home vs. Counties Without") +
  xlab("Days Since Stay-at-Home") +
  coord_cartesian(xlim=c(min(deaths$stay_day), max(deaths$stay_day)), ylim=c(min(deaths$lo), max(deaths$hi)), expand=TRUE) +
  theme(
        axis.title.x = element_text(size=20, angle=0, face="bold"),
        axis.title.y = element_text(size=21, angle=90, margin = margin(t = 0, r = 10, b = 0, l = 0), face="bold"),
        axis.text.y = element_text(size=13, angle=0, face="bold"),
        axis.text.x = element_text(size=13, face="bold"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line(),
        legend.title=element_blank(),
        legend.position = c(.2, 0.9),
        legend.text=element_text(size=14, face="bold"),
        legend.key = element_rect(colour = "black", fill="transparent")
)

# plot----
png(filename = 'figure4.png', width=15, height=8, units = "in", type="cairo", res=300)
grid.arrange(casesplot, deathsplot, ncol=2)
grid.text("a", x=unit(0.085, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
grid.text("b", x=unit(0.585, "npc"), y=unit(0.98, "npc"), gp=gpar(fontsize=30,font=2))
dev.off()
```
